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Abstract — Radio  frequency  (RF)  tomography  is  proposed  to  de¬ 
tect  underground  voids,  such  as  tunnels  or  caches,  over  relatively 
wide  areas  of  regard.  The  RF  tomography  approach  requires  a 
set  of  low-cost  transmitters  and  receivers  arbitrarily  deployed  on 
the  surface  of  the  ground  or  slightly  buried.  Using  the  principles 
of  inverse  scattering  and  diffraction  tomography,  a  simplified 
theory  for  below-ground  imaging  is  developed.  In  this  paper, 
the  principles  and  motivations  in  support  of  RF  tomography 
are  introduced.  Furthermore,  several  inversion  schemes  based 
on  arbitrarily  deployed  sensors  are  devised.  Then,  limitations  to 
performance  and  system  considerations  are  discussed.  Finally,  the 
effectiveness  of  RF  tomography  is  demonstrated  by  presenting 
images  reconstructed  via  the  processing  of  synthetic  data. 

Index  Terms — Buried-object  detection,  ground-penetrating 
radar  (GPR),  inverse  scattering,  radio  frequency  (RF)  tomogra¬ 
phy,  tunnel  detection. 

I.  Introduction 

IN  RECENT  years,  underground  void  detection  and  local¬ 
ization  has  become  a  critical  task,  motivated  by  the  incum¬ 
bent  need  for  protecting  national  borders  and  for  monitoring 
sensitive  areas,  such  as  prisons,  banks,  and  power  plants.  Ad¬ 
ditionally,  tunnel  detection  is  imperative  for  civil  applications, 
including  mining  safety,  search  and  rescue  in  devastated  ar¬ 
eas,  environmental  engineering,  geophysics,  archaeology,  and 
speleology. 

Presently,  information  concerning  voids  beneath  the  ground 
is  typically  obtained  using  different  techniques,  such  as  micro¬ 
gravity  (MG)  [1]— [3],  electrical  resistivity  tomography  (ERT) 
[4],  [5],  seismic  sensing  [6]-[8],  magnetotellurics  (MT)  [2], 
[9],  and  ground-penetrating  radars  (GPRs)  [10],  [11].  Unfor¬ 
tunately,  MG  is  suited  for  shallow  targets,  while  MT  is  more 
appropriate  for  very  deep  targets  [9] ;  ERT  may  not  work  ade¬ 
quately  when  the  impedance  of  the  ground  is  high,  and  seismic 
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methods  suffer  from  strong  external  noise,  lack  of  polarization 
diversity,  weathered  soils,  and  high  attenuation  [12]. 

GPR  appears  as  the  most  versatile  approach.  However,  it 
cannot  yet  be  considered  a  reliable  and  practical  method  in 
tunnel  detection  for  the  following  reasons.  First,  the  resolution 
for  classical  GPRs  is  generally  improved  by  using  large  band¬ 
width,  often  requiring  high  frequencies.  However,  due  to  soil 
properties,  higher  frequencies  experience  higher  attenuations, 
and  the  increased  frequency/bandwidth  affects  the  signal-to- 
noise  ratio  (SNR)  and  intensifies  dispersion  effects.  Second, 
systems  using  lower  frequencies  [13],  [14]  require  electrically 
small  wideband  antennas,  which  results  in  complex  wideband 
systems;  yet,  they  may  not  provide  adequate  resolution  to  de¬ 
termine  the  geometry  of  targets.  Third,  the  available  frequency 
spectrum  for  some  applications  can  be  severely  limited  by 
interference  with  external  sources  (e.g.,  broadcasting  stations); 
therefore,  a  reliable  system  must  operate  using  a  small,  discrete, 
and  selectable  number  of  frequencies.  Fourth,  the  interpretation 
of  the  raw  data  is  affected  by  the  operator’s  expertise,  and 
a  priori  information  is  necessary  to  obtain  reliable  results  [10]. 

To  overcome  attenuation  losses  (thus  enabling  detection  at 
large  investigation  depths),  GPR  may  be  inserted  into  boreholes 
[12]-[24],  which  generally  provides  an  image  of  a  vertical 
section  in  the  plane  between  the  logs.  However,  for  the  impor¬ 
tant  task  of  tracing  tunnel  pathways  and  localizing  adits,  the 
“horizontal”  prospecting  is  more  desirable,  leaving  the  deter¬ 
mination  of  depth  as  second  priority.  Furthermore,  boreholes 
are  expensive,  subject  to  drilling  misalignments,  and,  most 
important,  impractical  in  inaccessible  terrains. 

The  enhancement  of  resolution  may  be  accomplished  by  im¬ 
plementing  the  principles  of  RF/microwave  tomography  in  the 
framework  of  GPR  [25]-[31],  referred  to  as  GPRT.  Although 
an  improvement  in  reconstructed  images  may  be  achieved, 
GPRT  still  fails  to  be  reliable  for  tunnel  detection  for  several 
reasons,  including  as  follows:  1)  GPRT  is  typically  employed 
in  a  bistatic  configuration,  where  receive  and  transmit  antennas 
are  separated  by  an  electrically  small  distance,  while  operating 
in  proximity  to  the  air-earth  interface;  therefore,  it  suffers 
from  limited  view  diversity  (using  finite  observation  domains) 

[32] .  2)  To  compensate  for  this  limited  view,  an  increase  in 
the  amount  of  information  is  achieved  by  using  multitude  of 
frequency  tones,  thus  making  the  system  wideband  [28],  [31], 

[33] .  3)  Any  GPR  survey  requires  the  operator  to  be  above  the 
area  under  investigation.  To  date,  no  system  has  been  designed 
to  work  effectively  in  remote  situations  where  one  is  disallowed 
to  reach  the  area  of  regard. 

This  paper  extends,  further  improves  upon,  and  gives  a 
unified  framework  to  preliminary  ideas  in  RF  tomography 
proposed  in  [33],  [37]-[42].  In  particular,  Section  II  describes 
the  concept  of  RF  tomography  for  belowground  imaging; 
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Fig.  1 .  Principles  of  RF  tomography  for  tunnel  detection.  Transmitters  radiate 
power  into  the  ground.  Receivers  collect  the  scattered  field  and  send  this 
information  to  the  main  station. 

Section  III  illustrates  the  principles  of  RF  tomography, 
presenting  the  related  forward  model.  Section  IV  describes 
several  inversion  schemes.  Finally,  Section  V  shows  several 
reconstructed  images  obtained  via  the  inversion  techniques 
described  in  Section  IV. 

II.  RF  Tomography 

To  address  the  unresolved  problems  in  Section  I,  we 
introduce  a  new  approach  based  upon  a  multiview/multistatic 
design,  where  the  view  (associated  with  the  distributed 
transmitters)  and  the  observation  (associated  with  multiple 
receivers)  diversities  increase  the  information  pertaining  to 
the  scene  [34],  [35]  while  reducing  the  necessity  for  a  large 
spectral  content  (bandwidth)  of  the  probing  waveform:  In 
principle,  with  just  a  single  frequency,  it  is  possible  to  obtain 
high-resolution  images  [34],  [36]. 

The  proposed  approach  considers  two  separate  sets  of  N 
electromagnetic  transmitters  (Tx)  and  M  electromagnetic  re¬ 
ceivers  (Rx),  commonly  referred  with  the  generic  name  of 
Transponder.  These  transponders  are  placed  on  (or  in)  the 
ground  at  an  arbitrary  position  properly  defined  by  the  oper¬ 
ator.  The  transmitter  Tx  radiates  a  known  waveform  using  a 
suitable  polarization.  The  probing  wave  impinges  upon  a  buried 
dielectric  or  conductive  anomaly,  thus  generating  a  scattered 
wave-field.  The  receivers  collect  samples  of  the  scattered 
instantaneous  electric  field  and  estimate  the  complex-valued 
electric  field  phasor  at  their  locations.  Subsequently,  this  in¬ 
formation  is  relayed  to  a  data  collector  (Fig.  1).  To  ease  the 
detection  process,  one  transmitter  and  frequency  of  operation 
per  time  are  activated,  thus  simplifying  the  receiver’s  capability 
to  properly  discern  the  origin  of  the  incoming  wave-field.  For  a 
given  sampling  time,  the  used  spectrum  is  virtually  restricted  to 
a  single  frequency  component,  thus  ensuring  ultranarrowband 
system  architecture,  low  noise,  and  affordable  cost.  To  ease  the 
set  up  and  portability  of  the  system,  sensors  are  intended  to  be 
“dart”  shaped,  as  shown  in  Fig.  2.  Sensors  may  be  equipped 
with  built-in  GPS  for  precision  timing  and  positioning  and  an 
S-band  communication  link  to  transfer  the  collected  data  to  the 
overhead  base  station.  During  the  sensor-dropout  process,  some 
transponders  may  fall  in  obstructed  or  heavily  cluttered  regions 
(e.g.,  Tx  and  Rx  lay  in  proximity  or  over  a  vegetation  layer)  or 
they  may  not  survive  the  impact  with  the  ground.  In  any  case, 
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Fig.  2.  Proposed  system  design  for  transponders.  The  “dart”  shape  facilitates 
the  penetration  in  the  ground.  Not  to  scale. 

the  proposed  reconstruction  process  accounts  for  the  eventual 
failure  or  obstruction  of  transponders,  by  properly  neglecting 
corrupted  sensors. 

This  approach  may  also  operate  using  a  discrete  set  of  mono¬ 
chromatic  frequency  components,  selected  according  to  envi¬ 
ronmental  conditions.  A  suitable  modulation  is  stepped  FM; 
however,  in  this  paper,  the  conclusions  are  derived  indepen¬ 
dently  from  the  type  of  modulation.  The  operating  frequencies 
must  allow  the  electromagnetic  wave  to  penetrate  deeply  into 
the  ground,  while  simultaneously  provide  acceptable  detection 
and  resolution  capabilities.  Higher  frequencies  lead  to  better 
resolution,  but  strong  attenuation  limits  the  range  of  operation 
[60] .  Conversely,  if  the  frequencies  are  low,  the  corresponding 
resolution  may  not  be  adequate  to  localize  tunnels,  and  the  field 
behavior  becomes  diffusive,  thus  reducing  the  backscattered 
field  [29].  Based  upon  electrical  parameters  of  rocks  reported 
in  [10],  [61],  and  [62],  a  suitable  range  of  frequencies  for  this 
application  is  the  range  of  1-15  MHz,  but  the  final  choice 
strictly  depends  on  the  expected  target  type  and/or  depth. 

When  antennas  are  located  belowground,  three  different 
modes  of  propagation  between  transmitter  and  scatterer,  and 
between  scatterer  and  receiver,  are  excited:  direct,  reflected,  and 
lateral  waves  [52]-[54].  Generally,  the  lateral  mode  of  propaga¬ 
tion  is  the  most  undesirable:  although  it  could  be  estimated  and 
removed,  its  contribution  to  the  overall  field  may  be  so  high  that 
it  may  saturate  the  LNA  at  the  receiver  side,  thus  masking  the 
weak  signal  coming  from  the  scatterer. 

Using  vertical  dipoles  [40],  [56],  vertical  ferrite-loaded  coils 
[55],  or  gradiometer  antennas  [18],  [41],  [52],  [55],  [58],  the 
effect  of  lateral  wave  may  be  reduced  [42].  However,  recall 
that  the  proposed  method  is  valid  in  principle  for  any  envi¬ 
ronment,  and  the  reconstruction  accuracy  depends  primarily  on 
the  proper  choice  of  the  Green’s  function  (see  next  section  for 
details). 


III.  Forward  Model 

The  first  step  to  obtain  a  general  tomographic  inversion 
procedure  is  to  define  a  suitable  model  for  the  electric  field  that 


Authorized  licensed  use  limited  to:  WASHINGTON  UNIVERSITY  LIBRARIES.  Downloaded  on  August  16,2010  at  23:16:23  UTC  from  IEEE  Xplore.  Restrictions  apply. 


1130 


IEEE  TRANSACTIONS  ON  GEOSCIENCE  AND  REMOTE  SENSING,  VOL.  48,  NO.  3,  MARCH  2010 


Fig.  3.  Three-dimensional  geometry  for  the  model. 

closely  represents  the  actual  scene  while  keeping  enough  sim¬ 
plicity  to  facilitate  the  subsequent  inversion.  Similar  procedures 
are  also  discussed  in  [2],  [44],  and  [46]. 

The  3-D  geometry  shown  in  Fig.  3  is  considered.  For  sim¬ 
plicity,  a  single  operating  frequency  /  is  assumed,  but  extension 
to  the  multifrequency  operation  is  straightforward  [40].  Under 
the  monochromatic  assumption,  the  ground  is  modeled  as  a 
homogeneous  medium  with  relative  dielectric  permittivity  £p, 
conductivity  op,  and  permeability  fi o.  The  targets  (i.e.,  tunnels 
or  voids)  are  assumed  to  reside  in  the  investigation  domain  D. 
The  sources  are  N  electrically  small  dipoles  (of  length  Alf)  or 
loops  (of  area  A*)  fed  with  current  Jt  and  located  at  position 
(view  diversity)  and  with  (electric  or  magnetic)  dipole  moment 
directed  along  the  unit  vector  a^.  For  each  transmitting  antenna, 
the  scattered  field  E s  is  collected  by  M  receivers  (observation 
diversity),  located  at  points  in  space. 

We  assume  the  relative  dielectric  permittivity  profile  er(rf) 
and  the  conductivity  profile  cr(r')  inside  the  investigation  do¬ 
main  D  as  unknowns  of  the  problem. 

Accordingly,  the  inverse  problem  is  recast  in  terms  of  the 
unknown  dielectric  permittivity  contrast  [35],  [43] 

e5(r')  =  er(r')  -  £d  +  j'A-!-  (1) 

27r/£o 

In  this  way,  the  wavenumber  inside  D  can  be  expressed  as  [30] 

k2(  r')  =u2no£o£r(r')  +joJH0a(r') 

=  k2D  +  klss(r')  (2) 

kD  =  wyV0£ 0£d 

k0=u£iM^.  (3) 

The  function  in  (1)  accounts  for  the  difference  between  the 
unknown  dielectric  permittivity  of  the  object  and  that  of  the 
host  medium.  We  also  accounted  for  the  conductivity  profile 
in  (1)  and  (3),  because,  in  general,  the  halo  of  a  tunnel  shows 
an  increased  apparent  conductivity  [12],  primarily  due  to  salt 
dissolved  in  wet  surfaces,  steel-reinforced  concrete  for  ground 
control,  power  lines,  rails,  ventilation  tubes,  etc. 

For  each  point  r'  in  region  D ,  the  vector  wave  equation  holds 

V  x  V  x  E(r')  =  [kl  +  k2es(r')]  E(r').  (4) 


The  scattered  wave  in  a  point  r  ^  D  that  is  a  solution  of  (4)  can 
be  written  in  terms  of  integral  equation  of  the  dyadic  Green’s 
function 


Es(r)  =  kl  JJJ  g±(r,r')  ■  E(r')es(r')dr'  (5) 

D 

where  E(r')  is  the  total  field  in  the  investigation  domain  D , 
given  as  the  superposition  of  the  incident  field  E1  ( r')  (i.e.,  the 
field  in  the  investigated  area  when  objects  are  absent)  and  the 
field  E5(r)  scattered  by  the  targets. 

The  inverse  scattering  problem  in  (5),  i.e.,  finding  £<5(1*') 
based  on  the  knowledge  of  Es(r),  is  nonlinear.  Nevertheless, 
it  can  be  recast  as  a  linear  problem  by  means  of  the  Bom 
approximation  (BA).  In  fact,  from  the  Lippman-Sch winger 
representation  [57]  of  the  scattering  problem,  the  total  field  in 
region  D  can  be  expressed  in  operator  form 

E(r)  —  E7(r)  =r^E(r) 

E(r)  =[I-r^]-1E/(r)  (6) 

where  I  is  the  identity  vector,  T  represents  the  operation  of 
convolution,  and  4/  represents  the  operation  of  multiplication 
by  eg.  If  the  operator  norm  [71]  of  is  less  than  one,  then  we 
can  expand  (6)  using  Neumann  series 

E(r)  =  (I  +  r^  +  r^r^  +  ---)E7(r).  (7) 

The  first-order  BA  accounts  only  for  the  first  term  in  (7),  thus 
the  total  field  inside  the  integrand  of  (5)  can  be  approximated 
by  the  known  incident  field  [64],  yielding 

Es(r)  =  kl  jjj g(r,  r')  •  E1  (rr)ss(r')dr' .  (8) 

D 

Therefore,  the  corresponding  inverse  problem  is  cast  as  the 
inversion  of  the  linear  integral  equation  connecting  the  permit¬ 
tivity  contrast  function  to  the  scattered  field  data. 

The  use  of  BA  is  justified  by  considering  the  following 
reasons. 

1)  Tunnels  (or  other  targets  of  interest)  are  isolated,  limited 
in  number,  and  embedded  in  a  lossy  medium.  Therefore, 
mutual  interaction  between  tunnels,  a  phenomenon  ne¬ 
glected  by  BA,  can  be  assumed  negligible. 

2)  In  general,  the  inhomogeneities  of  the  soil  are  electrically 
small,  and  their  conductivity  remains  low.  Therefore, 
their  scattered  fields  are  insignificant  as  compared  to  the 
RF  signal  reirradiated  by  tunnels. 

3)  The  operator  norm  of  T4/  is  generally  <  1,  due  to  the 
losses  in  the  ground  (that  introduce  an  exponential  atten¬ 
uation  of  the  fields)  [68]  and  to  the  spherical  spreading 
factor  of  the  incident  and  scattered  waves. 

4)  The  main  goal  is  to  detect,  localize,  and  approximately 
determine  the  geometry  of  the  targets.  Toward  this  ob¬ 
jective,  BA-based  inversion  algorithms  preserve  the  in¬ 
formation  on  target  localization,  even  when  objects  are 
strong  scatterers  (although  the  quantitative  description  of 
es  in  D  is  generally  prejudiced)  [30],  [43]. 
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The  incident  field  can  be  expressed  in  terms  of  Green’s 
functions 


E«  (A  A)  =  QQ  (A  A)  •  A  (9) 

where  Q  =  juj^AFP  for  an  electrically  small  dipole  or  Q  = 
—juj/ivA1!1  for  an  electrically  small  loop. 

Additionally,  the  field  received  by  a  dipole  or  loop  with 
moment  direction  positioned  at  due  to  an  equivalent  (in 
terms  of  E^)  current  distribution  defined  inside  the  investiga¬ 
tion  domain  D  can  be  expressed  as  [28] 

£S(A -Ai)=fc0  JJJ  Ai  -£Am>A-En  (A  A)  £s(r')dr'. 

D 

(10) 

Substituting  (9)  in  (10),  we  obtain  the  scalar  forward  model 
for  the  scattered  field  as 

es  (A,  A 0  =L  <A(A) 

=  Qk o  JJJ  [am  '  A  ( A>  r')] 

•  lii(A  A)  •  A]  ^(r')dr'.  (11) 

Equation  (11)  shows  a  linear  relation  between  scattered  field 
and  the  dielectric  profile,  and  the  reconstruction  procedures 
presented  in  Section  IV  are  based  on  its  inversion. 

However,  the  actual  electric  field  measured  at  the  receiver 
is  the  superposition  of  the  scattered  field  from  targets  and  the 
direct  coupling  between  Tx  and  Rx.  From  (11)  and  (9),  we 
determine  the  total  field  experienced  at  the  receiver  in  time 
domain 


+A,A+ 

=  Re  I  e~jut 


Qko  JJJ  [Ai-£(Ai,A] 

•  [£  (A  A)  •  A]  £s (Adr'  +  QA 


•£(rm-A)  •  A  +  #  AAA) 


(12) 


The  term  •  G(r^,  rJJ  •  in  (12)  represents  the  direct 
coupling  between  Tx  and  Rx,  and  it  can  be  considered  a  source 
of  deterministic  clutter.  The  term  iT(rJ^,r^)  represents  the 
nonlinear  contribution  due  to  higher  order  Born  series  terms. 
In  the  reconstruction  process,  H( rJJ  is  considered  as  un¬ 
known  bias  (clutter)  to  the  desired  signal.  The  random  variable 
N(t)  can  be  modeled  as  Gaussian  process  with  zero  mean  and 
variance  equal  to  the  noise  power 


PN  =  Fa  +  10  log  5  +  10  log [KbT0]  [dB  •  W]  (13) 


where  B  is  the  bandwidth,  Kb  is  the  Boltzmann  constant,  T0 
is  the  environmental  temperature,  and  the  external  noise  figure 
Fa  can  be  inferred  by  consulting  [63]. 

Tomography  is  inherently  suited  for  noise  mitigation,  since 
it  is  ultranarrowband  (therefore,  Pn  is  intrinsically  very  low), 
and  by  simply  averaging  n  samples  of  the  same  signal,  we 


obtain  a  theoretical  SNR  increase  of  yjn.  Furthermore,  the  static 
contributions  in  (12)  as  functions  of  r^,  are  generally  low 
correlated  with  the  value  of  the  scattered  field,  meaning  that  the 
combination  of  view  and  observation  diversities  randomizes  the 
clutter  due  to  the  direct  path  coupling. 

To  obtain  the  electric  field  in  phasor  form  from  (12),  the 
received  field  can  be  mixed  with  two  coherent  oscillators  to 
retrieve  the  in-phase  and  quadrature  components,  since 

E  (A-  Ax,  *)  =  Re  [E  (A,  A)  exp  (-jut)] 

—  Re  [E  (A,  A)]  cos (ut) 

+ Im  [E  (A.  A)]  sin (14) 

Hence,  the  real  and  imaginary  parts  are 

t+T 

Re  [E  (A,  rrn)]  =J  J  E  (Ai  A,  t)  cos (ut)dt  (15) 

t 

t+T 

Im  [E  (A,  A)]  =  |;  J  E  (A,  A,  t)  sin (wt)dt.  (16) 

t 

IV.  Inversion  Procedures 

From  a  mathematical  point  of  view,  the  problem  of  finding 
the  dielectric  profile  is  to  compute  the  inverse  of  the  linear  op¬ 
erator  L  in  (11),  connecting  the  unknown  dielectric  profile  and 
the  scattered  field  data.  Noise  and  static  clutter  contributions 
are  assumed  sufficiently  small,  so  that  they  can  be  considered 
as  perturbations  on  the  measured  data. 

A  way  to  compute  L-1  is  to  perform  a  numerical  inversion 
of  L.  Fet  us  collect  the  sampled  field  data  in  an  ordered  NM 
vector  E_s  =  {Es( r^,rj^)}  and  discretize  the  domain  region 
D  in  K  voxels,  each  one  located  at  position  r'k:  The  contrast 
dielectric  permittivity  can  be  embodied  in  a  column  vector  e  §  = 
{^(r^)}  of  length  K ,  and  it  represents  the  set  of  unknown 
parameters.  After  this  discretization,  (11)  can  be  rewritten  in 
matrix  form 

Es  =  h£5  (17) 

where  L  is  now  a  matrix  with  dimensions  NM  x  K  and  Es, 
e§  are  column  vectors. 

The  problem  is  then  to  invert  the  relation  (17).  Since  L  is 
generally  not  a  square  matrix,  we  need  to  consider  its  pseudoin¬ 
verse,  which  we  will  still  indicate  with  L-1.  Due  to  the  usual 
location  of  Tx,  Rx,  and  targets,  L  is  generally  ill  conditioned.  A 
common  way  to  quantify  the  behavior  of  L  is  by  inspection  of 
its  condition  number  k.  For  the  operator  L,  it  is  quite  common 
to  obtain  values  of  k  above  106. 

This  leads  to  artifacts  in  the  reconstruction  process,  particu¬ 
larly  exacerbated  when  noise  (thermal,  external,  quantization) 
or  clutter  is  impinging  upon  receivers. 

According  to  the  accuracy  required  from  the  system,  we 
present  four  inversion  strategies. 

1)  Levenberg-Marquardt  (FM)  regularization  procedure. 
This  method  is  relatively  accurate  for  any  environmental 
condition,  and  it  is  robust  in  presence  of  noise.  It  requires 
a  proper  choice  of  a  regularization  parameter. 
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2)  Truncated  singular-value  decomposition  (TSVD).  This 
method  is  also  relatively  accurate  in  any  scenario  and 
fairly  resistant  to  noise  interference.  This  method  also  of¬ 
fers  deeper  insight  into  the  physics  behind  the  reconstruc¬ 
tion,  and  the  output  can  be  easily  adjusted  by  properly 
selecting  the  number  of  meaningful  singular  values.  In 
addition,  the  number  of  retained  singular  values  in  TSVD 
plays  the  same  role  of  the  LM  regularization  parameter. 

3)  Back-propagation  approach.  This  method  works  properly 
only  when  the  operator  L  is  relatively  well  conditioned. 
This  implies  that  it  can  be  used  only  for  particular  config¬ 
urations  and  when  the  SNR  is  relatively  high.  However, 
the  computational  time  is  drastically  reduced. 

4)  Fourier-Bojarski  approach.  This  is  the  fastest  inversion 
scheme,  and  it  is  suited  for  far-held  probing  and  quasi¬ 
lossless  Green’s  function.  It  privileges  speed  instead  of 
accuracy. 

A.  LM  Regularization 

An  efficient  method  to  compute  the  inverse  of  an  ill- 
conditioned  matrix  is  by  using  the  LM  regularization  procedure 
[2].  In  this  way,  the  dielectric-permittivity  contrast  is  esti¬ 
mated  as 

!*(/?)  =  (LHL  +  /3I  )~%HES  (18) 

where  LH  denotes  the  adjoint  of  L  and  (3  is  the  regular- 
ization  parameter  in  the  Tikhonov  sense,  which  needs  to  be 
advantageously  selected.  Since  a  proper  choice  of  /3  may  be  a 
difficult  task,  our  initial  guess  for  [3  is  the  midpoint  value  of  the 
singular- value  dynamic  range  of  L.  Oftentimes,  it  is  necessary 
to  determine  /3  through  optimizations  and  iterations  before  a 
meaningful,  sharp,  and  low-blurred  image  is  reconstructed  [2], 
[69].  This  implies  that  a  (computationally  expensive)  matrix 
inversion  for  each  attempt  may  be  necessary.  To  accelerate 
this  process,  the  SVD  may  be  used  advantageously,  which  is 
described  next. 

B.  Truncated  SVD 

A  more  efficient  way  to  invert  the  ill-conditioned  matrix  L  is 
proposed  in  [34],  [35],  and  [43]  and  takes  advantage  of  SVD. 
In  fact,  L  can  be  decomposed  as  follows: 

L  =  H§YH  (19) 

where  S  is  a  diagonal  matrix  containing  the  ordered  singular 
values  Si  of  L.  The  pseudoinverse  of  L  can  be  written  as 

If1  =  Yi_1UH.  (20) 

Singular  values  associated  to  S-1  that  are  considerably  large 
as  compared  to  1/si  represents  the  sensitive  directions  of  L: 
along  these  directions,  small  amounts  of  noise  or  clutter  in  the 
collected  electric  field  will  lead  to  a  large  (undesired)  devi¬ 
ation  of  £§. 

A  way  to  remove  this  sensitivity  is  to  consider  only  the  first  k 
smaller  singular  values  of  S_1  and  setting  to  zero  the  remaining 
large  ones.  This  strategy  is  commonly  referred  to  as  'Truncated 


SVD ”  [66].  Accordingly,  the  dielectric  profile  can  be  esti¬ 
mated  as 

ls(k)=Yg~1  nH  (21) 

The  SVD  can  also  be  very  useful  to  properly  dimension  the  [3 
parameter  in  the  LM  method.  In  fact,  if  we  rewrite  (18)  in  terms 
of  (19),  we  obtain 

If1  =  (LHL  +  /5I)_1Lh  =  Ydiag  (^L-2)  (22) 

\si  +  [3  Si  J 

TSVD  and  the  SVD  representation  of  the  LM  method  in  (22) 
have  a  remarkable  feature.  In  both  cases,  once  the  evaluation  of 
the  singular  system  of  L  is  performed  (19),  the  reconstructed 
image  is  simply  computed  by  a  (fast)  matrix  multiplication  as 
in  (20).  This  means  that  a  new  image  is  obtained  by  varying  the 
number  of  singular  values  (for  the  TSVD  method)  or  varying 
(3  (for  the  LM  method),  and  it  can  be  computed  extremely  fast 
(often  in  real  time). 

C.  Back  Propagation 

The  approximate  number  of  computations  to  perform  the 
SVD  of  L  is  9 K3  +  12 MNK2  [47].  When  the  size  of  the 
matrix  is  large,  the  evaluation  of  (19)  may  become  prohibitive, 
and  an  alternative  efficient  strategy  has  to  be  pursued. 

In  this  case,  the  contrast  permittivity  function  can  be  esti¬ 
mated  using  the  following  approximation: 

Is  =  (LHL)_1LH=ES  =  LH£S.  (23) 

Therefore,  the  inversion  is  carried  out  by  just  using  the  adjoint 
of  L.  This  particular  inversion  holds  when 

LHL  =  a  I  a  >  0  (24) 

where  al  represents  a  scaled  identity  matrix.  This  result  implies 
that,  in  principle,  ft(LHL)  =  1.  However  [47] 

k(LHL)  =  4LH)2  =  4L)2.  (25) 

Therefore,  the  use  of  the  adjoint  for  the  inversion  is  acceptable 
only  when  the  singular  values  of  L  exhibit  a  limited  dynamic 
range.  In  this  case,  an  explicit  formula  for  the  solution  of  (23) 
in  terms  of  each  r'k  can  be  formulated 

N  M 

£<5  (4)  =  Q*k o  xm  (a«  '  =H  (4».  4)) 

n=  1  m=  1 

•  (4> 4)  •  4)  Es  (4,  rrm ) .  (26) 

Equation  (26)  is  commonly  referred  to  as  matched  filtering , 
migration ,  or  hack  propagation  [67].  This  technique  is  suited 
for  parallel  processing  [33]. 

D.  Fourier-Bojarski  Approach 

A  simple  and  fast  approach  (albeit  less  accurate)  is  to  take 
advantage  of  the  Fourier  relation  arising  between  scattered  field 
and  object  shape,  as  discussed  in  the  literature  under  the  topic 
of  diffraction  tomography. 
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In  fact,  if  targets  and  sensors  are  distant  enough  so  that  the 
EM  field  can  be  approximated  as  a  locally  lossless  plane  wave 
(normally  occurring  when  the  fields  are  primarily  propagating 
as  1/r),  then  the  forward  model  can  be  expressed  as  explained 
as  follows. 

We  define  the  unit  norm  direction  of  propagation  vectors  as 


I*n  =  —  x  sin  Oj  cos  (fn  —  y  sin  0 J  sin  y>ln  -  z  cos  0 ln  (27) 

lm  =  x  sin  6rm  cos  tprm  +  y  sin  9rm  sin  <pTm  +  z  cos  6>4 .  (28) 

Using  the  paraxial  approximation,  the  transmit  Green’s  func¬ 
tion  at  the  generic  position  r'  inside  region  D  can  be  sim¬ 
plified  as 


Gn(  <y)  = 


exp 


(+jkDrtn)  exp  (+jkD Pn  ■  r ') 
47r  ri 


(29) 


while  the  receive  Green’s  function  can  be  expressed  as 

exp  ( +jkDrrm )  exp  [~jkD Ym  •  i-7) 


Grm(  r'ry  = 


47r  r!! 


(30) 


Therefore,  for  a  pair  of  transmitters  and  receivers,  the  scattered 
field  can  be  rewritten  as  [36] 


Es  (r1  rr  )  - 
V  n’  m)  — 


1,2  z.t  .  Ar 
^0an  am 

IOtt2?^?! 


Qe+j&r>(rr+rm) 


III £s(' r  ) exp  0"  “  1™)  ■  ^  dr ’ '  (31) 


D 


The  quantity  kr>  (1^  —  1J^)  can  be  represented  by  a  3-D  vector 


kD  (if,  ~  t)  • 


(32) 


Equation  (31)  is  then  rewritten  as  [33] 


Es(kmn)  = 


^0dn  ’ 


167r2r^r!l 


ge+i/c0(r^+r;) 


JJJ £S{ r')  exp(+jkmra  •  r')dx' .  (33) 

D 


It  is  useful  to  consider  a  normalized  version  of  (33) 


_ § 

approximated  (in  the  limit)  as  a  continuous  function  E  (k), 
and  (34)  can  be  interpreted  as  a  3-D  inverse  Fourier  transform 
of  the  permittivity  contrast  function.  Therefore,  an  image  can 
be  reconstructed  via  direct  Fourier  transformation  of  (33),  i.e., 

£$(r')  =  JJJ  Es( k)  exp(— jk  •  r')dk  (35) 

K 

_ g 

where  the  domain  of  integration  K  is  the  support  of  E  (k). 
By  inspection  of  (32),  we  conclude  that  when  the  sensors 
completely  encircle  the  target,  K  is  a  sphere  of  radius  2k d-  In 
other  words,  the  available  information  of  the  spectral  content  of 
^(r')  is  limited  up  to  the  spectral  component  contained  inside 
a  sphere  of  radius  2kp-  Therefore,  the  reconstructed  image  of 
the  dielectric  profile  will  be  a  low-pass-filtered  version  of  the 
true  image. 

By  studying  the  impulse  response  of  (34),  the  minimum  res¬ 
olution  achievable  using  Fourier-Bojarski  approach  is  shown 
to  be  d  =  Ad/3  [39].  For  half-space  problems,  the  resolution  is 
further  reduced  [49]. 

In  real  cases,  where  a  finite  number  of  sensors  are  deployed 
(i.e.,  the  spectral  domain  is  undersampled)  and  external  noise 
affects  the  measurements,  the  resolution  is  lower,  and  artifacts 
in  the  reconstructed  image  are  very  common.  Severe  smearing 
and  blurring  effects  originates  mainly  from  the  invalidity  of  the 
paraxial  approximation.  A  way  to  overcome  this  limitation  is 
to  segment  the  region  D  in  smaller  analysis  regions  and  con¬ 
sider  an  inverse  problem  for  each  subregion,  and  the  resulting 
subimages  are  concatenated  to  form  the  final  image  [39].  The 
undersampling  of  the  spectral  domain  can  be  corrected  with 
several  approaches,  such  as  trilinear  interpolation  [50]  of  the 
available  samples,  or  using  projection  on  convex  sets  [51]  and 
smoothing  of  the  peaks  in  the  spectral  domain  to  estimate  the 
missing  samples  [39]. 


V.  Simulations 

We  present  some  simulation  results  to  test  the  proposed 
belowground  imaging  system.  A  set  of  six  transmitters  and 
26  receivers,  operating  at  5  MHz,  is  placed  belowground  at 
a  depth  d  =  0.15  m  (a^  =  ar  =  x),  as  shown  in  Fig.  4.  Two 
empty  cylindrical  structures  with  radius  p  =  1  m  (representing 
two  tunnels)  are  assumed  to  be  embedded  in  a  host  medium 
with  relative  dielectric  permittivity  £d  =  10  and  conductivity 
aD  =  5  x  10-4  S/m  [62]  at  a  depth  h  =  25  m.  The  correspond¬ 
ing  attenuation  can  be  computed  using  [65] 


ES  (k  )  ~  rnrm  -jknjrj+r^)  fa  \ 


Qklki  •  a£ 

JJJ  es(  r)  exp(+jkmn  •  v')dr' .  (34) 


This  result  can  be  interpreted  in  the  following  way:  Each  col- 

_ § 

lected  sample  E  (kmn)  returns  the  value  of  the  kmn  spectral 
component  of  the  contrast  function  £<5(1*').  Theoretically,  if 
enough  samples  are  available  to  fully  populate  the  spectral 

_ g 

representation  of  ^(r'),  the  discrete  function  E  (kmn)  can  be 


a  =  4.34cr_Dv//io/(eo£D)  [dB/m],  (36) 


The  scattered  field  was  synthesized  using  the  FDTD  simu¬ 
lator  GPRMAX  [70]:  The  instantaneous  scattered  electric  field 
has  been  correlated  using  (15)  and  (16)  to  retrieve  the  electric 
field  in  phasor  form.  For  the  sake  of  speed  and  simplicity, 
in  (11),  we  selected  the  Green’s  dyad  for  the  homogeneous 
medium  with  the  same  properties  of  the  ground,  i.e.,  [44] 


£(r,r') 


i4'v 

K D 


ejkD  |r-r'| 

47r|r  —  r'| 


(37) 
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Fig.  4.  (Top  view)  Geometry  for  the  simulation. 


This  assumption  is  reasonable  when  the  sensors  are  deployed 
below  the  air/ground  interface,  the  frequencies  involved  are 
relatively  low,  and  targets  are  assumed  to  reside  nearly  perpen¬ 
dicular  to  sensors.  The  deeper  the  scatterer,  the  lower  the  lateral 
mode.  Nonetheless,  more  accurate  models  can  take  into  account 
the  distortions  due  to  half- space  or  layered  geometry  by  simply 
selecting  the  proper  Green’s  function  under  a  spectral  form 
[28],  [29],  [44]-[46]  or  by  numerically  computing  the  Green’s 
function  using  method  of  moments  or  fast  eikonal  equation 
solvers  (subject  to  future  research). 

Direct-path  coupling  was  mitigated  thanks  to  the  exact 
knowledge  of  the  dyadic  Green’s  function  in  this  problem, 
which  enabled  accurate  determination  and  subsequent  cancel¬ 
lation.  Random  Gaussian  noise  has  been  added  on  the  data 
according  to  [63].  Nevertheless,  the  SNR  can  be  reduced  to 
a  desired  value  by  opportunely  sampling  and  averaging  the 
received  field,  as  discussed  in  Section  IV. 

Upon  completion  of  the  3-D  tomographic  inversion  (using  a 
mesh  of  1  m3  per  voxel),  the  horizontal  section  at  25-m  depth 
(constant  depth  slice)  is  plotted.  For  more  complex  scenarios,  a 
full  3-D  image  might  be  necessary. 

Dielectric  perturbations  of  values  less  than  5%  of  the  peak 
value  in  the  domain  D  have  been  neglected  in  the  final  images. 

The  inversion  procedure  using  Fourier-Bojarski  method  has 
a  different  (coarse)  voxel  size  equal  to  A^/4,  because  sampling 
at  finer  discretization  does  not  provide  more  information  (due 
to  resolution  limitations  [39]),  while  it  dramatically  increases 
the  sparsity  of  the  Fourier  domain,  leading  to  severe  artifacts. 

Fig.  5  shows  the  reconstructed  image  using  the  (fast) 
Fourier-Bojarski  method:  Although  the  resolution  is  coarse, 
basic  features  of  the  two  tunnels  can  be  discerned.  However,  the 
back-propagation  method  (see  Fig.  6)  is  clearly  showing  higher 
resolution  while  keeping  the  computational  cost  to  a  mini¬ 
mum.  For  high-level  image  reconstruction,  regularized  methods 
are  paramount.  In  Fig.  7,  the  image  has  been  reconstructed 
using  LM  method  (in  its  SVD  variant),  and  the  regularization 
parameter  has  been  empirically  selected  in  order  to  achieve 
the  sharpest  solution.  In  Fig.  9,  an  image  reconstructed  using 
TSVD  is  shown:  We  assume  that  10%  of  the  total  singular 


Fig.  6.  Reconstructed  image  using  back  propagation. 


X  Range  [m] 

Fig.  7.  Reconstructed  image  using  LM  method.  (3  has  been  empirically 
selected. 


values  represent  the  sensitive  directions  of  L,  and  therefore, 
they  will  be  not  included  in  the  reconstruction  based  on  TSVD 
(see  Fig.  8).  This  threshold  has  been  chosen  heuristically,  and  it 
may  vary  according  to  the  geometry  and  the  SNR. 
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Fig.  8.  Singular  values  behavior  (in  decibels  and  normalized  to  the  first 
singular  value)  of  the  L  operator. 


Fig.  9.  Reconstructed  image  using  TSVD.  The  number  of  singular  values  used 
has  been  empirically  selected. 


Finally,  Fig.  10  shows  the  reconstruction  starting  from  the 
scattered  field  data  computed  using  the  forward  model  (11)  and 
the  Green’s  function  in  (37).  Clearly,  the  TSVD  provides  very 
similar  reconstruction  results  for  the  two  cases  of  FDTD  and 
Born  field  data.  This  result  supports  the  accuracy  of  BA  for 
belowground  imaging  of  void  targets. 


VI.  Conclusion 

We  proposed  a  practical  method  for  tunnel  detection  that 
does  not  require  boreholes,  is  easily  deployable,  and  covers 
relatively  wide  areas. 

We  applied  diffraction  tomography  and  inverse-scattering 
principles  to  our  geometry.  We  proposed  four  simple  inversion 
methods  to  reconstruct  images  that  are  suited  for  the  environ¬ 
ment  encountered  in  practical  situations. 

Finally,  we  showed  several  reconstructed  images,  and  in  all 
cases,  the  location  of  the  two  tunnels  is  discernible  when  noise 
is  low,  while  tunnel  detection  amid  high  noise  is  only  possible 
when  LM  and  TSVD  methods  are  used. 

In  particular,  by  comparing  Figs.  6  and  7  with  Figs.  9  and 
10,  it  is  clear  that  LM  and  TSVD  methods  yield  higher  quality 
images.  Therefore,  if  the  computational  load  is  not  a  consid¬ 
eration,  these  two  methods  should  be  preferred  and  become 
imperative  when  environmental  conditions  are  hostile. 


Fig.  10.  Reconstructed  image  using  TSVD  and  scattered  field  data  generated 
using  the  forward  model  in  (1 1). 


The  proposed  strategy  offers  the  following  advantages. 

1)  RF  tomography  is  able  to  surveil  from  local/shallow 
to  global/deep  areas  of  regard,  and  rapidly  focus  on 
specified  areas,  by  simply  changing  frequency  of  oper¬ 
ation  and  the  delimitation  of  the  investigation  domain  D. 

2)  The  system  is  suited  for  both  cooperative  and  denied 
scenes,  where  the  physical  presence  of  the  human  oper¬ 
ator  is  hazardous. 

3)  The  ultranarrowband  sensor  design  and  the  system  archi¬ 
tecture  are  fiscal  and  manpower  affordable  solutions. 

4)  The  sensor  deployment  is  arbitrary  and  modular  (i.e.,  the 
addition  or  removal  of  sensors  does  not  compromise  the 
remainder  of  the  system). 

5)  The  resolution  of  the  system  can  be  sub  wavelength  and 
range  independent. 

However,  many  aspects  still  need  deeper  investigations,  such 
as  more  accurate  inversion  models,  the  use  of  different  Green’s 
functions,  improved  methods  of  direct-path  cancellation,  or 
more  considerations  on  the  actual  soil  and  antennas  behavior. 

In  particular,  the  issue  of  unwanted  lateral  waves  must  be 
addressed  in  two  ways:  by  investigating  clutter- suppression 
techniques  or  by  defining  a  suitable  Green’s  function  that  ac¬ 
counts  for  this  effect.  We  are  currently  pursuing  further  research 
in  both  directions. 
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